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ON A NON-REFLECTING BOUNDARY CONDITION FOR HYPERBOLIC 

CONSERVATION LAWS 

Ching Y. Loh 

Taitech, Inc. 

Beaver Creek, Ohio 45430 
Email: fsloh@turbot.grc.nasa.gov 

Abstract 

A non-reflecting boundary condition (NRBC) for practical computations in fluid dynamics and aeroacoustics 
is presented. The technique is based on the first principle of non-reflecting, plane wave propagation [1] and 
the hyperbolicity of the Euler equation system. The NRBC is simple and effective, provided the numerical 
scheme maintains locally a C 1 continuous solution at the boundary. Several numerical examples in ID, 2D 
and 3D space are illustrated to demonstrate its robustness in practical computations. 

1 Introduction 

It is well-known that non-reflecting boundary conditions (NRBCs) play an important role in fluid flow and aeroa- 
coustics computations. The need for artificial boundary conditions arises when the domain of the problem is un- 
bounded and extends to infinity. In order to treat the problem numerically, a domain of finite size is required and 
artificial boundaries are imposed. At these artificial boundaries, NRBCs are sought for to minimize their influences 
on the flow. A spurious reflection resulting from an inappropriate numerical boundary condition will contaminate the 
flow field and may entirely spoil the flow computation. Research on NRBC is a challenging topic in engineering and 
applied mathematics. For decades, a vast number of papers on NRBC have been published, e.g. see [3-6], the review 
paper by Givoli [2], and the references cited there. 

In one-dimensional flow, at an artificial boundary, Enquist and Majda [3], and Hedstrom [7] proved that a bound- 
ary condition is non-reflecting is equivalent to saying that the characteristic variables corresponding to the incoming 
characteristic curve remain constant across the artificial boundary (see also Hirsch [8], p. 370). For multi-dimensional 
flow, this 1-D technique is combined with dimension splitting and applied in the practical NRBC treatment. Such a 
combined treatment has been the topic of many papers on the characteristics based NRBCs (see e.g. [4] and references 
in [2]). 
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Other alternative treatments for NRBC found in the literature include (see [2-6]): 

(i) in the far-held, using a predictable asymptotic analytical solution at the boundary ( the radiation boundary 
condition), 

(ii) diminishing the strength of the waves/disturbances before they reach the artificial boundary, and thus minimizing 
the reflecting effect. Usually, increased numerical damping is applied to a zone between the core domain and the 
artificial boundary (the buffer zone or sponge zone) to do the job. 

(iii) In the recently developed PML (perfectly matched layer) method, a specially designed equation system is 
imposed in the matching layer (or sponge zone) to guarantee the exponential decaying of the disturbances in the layer 
[5,6], 

In the present paper, a different but simple criterion is introduced to treat the NRBCs of the time-dependent hyper- 
bolic conservation laws of gas dynamics. The criterion is based on the first principle of non-reflecting and plane wave 
propagation [1] rather than the characteristics theory. 

As it turns out that the NRBCs used in the recent CE/SE finite volume schemes for flow and aeroacoustics com- 
putations (e.g. [16,17]) can be directly derived from this criterion, the present paper also serves to explain why these 
simple NRBCs works well with the CE/SE schemes. 

The paper is arranged as follows: In Section 2, based on the first principle of non-reflecting, the propagation of plane 
wave and the hyperbolicity of the Euler equation system the continuity criterion of NRBC is introduced and proved. 
The numerical treatment of the continuity criterion is presented in §3, an extrapolation-like NRBC (Type I) based on 
this criterion and the numerical procedure are described. In §4 the relation between the NRBC and the flux balance 
across the boundary surface is established, which leads to another type of NRBC (Type II). Discussions on the NRBCs 
are given in §5. In Section 6, several numerical examples for outflow NRBC in one and multi-dimensional space are 
presented. Numerical examples with Type II NRBC at the inflow and other artificial boundaries are demonstrated in 
Section 7. Application of buffer/sponge zones is illustrated in Section 8. At last, the paper is concluded with remarks 
in Section 9. 

As the time-dependent hyperbolic conservation laws of gas dynamics ( in dimensionless form) is always incorpo- 
rated in the NRBCs in the present paper, they are briefly described here for later use: 

U t + F x + Gy + H z = Q, (1) 

where x, y, z and t are the streamwise and transversal coordinates and time, respectively. The conservative flow 
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variable vector U and the flux vectors in the streamwise and radial directions, F, G, and H, are given by: 
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where it, v, w and p, p are respectively the velocity components, density and pressure, e = + 1/2 (u 2 +v 2 +w 2 ), 

and the enthalpy H = p/p + e with 7 = 1.4. The right hand side Q is the source term which may include the possible 
external forcing terms and/or viscous fluxes. 

By considering (x, y, z, t) as coordinates of a four-dimensional Euclidean space, E 4 , and using Gauss’s divergence 
theorem, it follows that Eq. (1) is equivalent to the following integral conservation laws: 


I m ds= / Q m dF, m = 1,2, 3, 4, 5, 


( 2 ) 


JS(V) Jv 

where 5(F) denotes the surface around a volume V in E 4 and I m = (F m , G m , H m , U m ) stands for the flux vectors, 
I m ds = I m • mis , n being the outgoing unit normal vector in E 4 , is the flux at the infinitesimal surface element ds. 


2 The continuity criterion of NR6C 

In numerical schemes, particularly in finite volume numerical approaches with hyperbolic conservation laws, grid 
nodes are often cell centers and the boundary faces are often formed by the boundary cell surfaces. No node lies exactly 
on the boundary. As such, a continuity criterion of NRBC and the consequent NRBC treatments are introduced based 
on the first principle of non-reflecting, the first principle of plane wave propagation [1] and the hyperbolicity of the 
equation system. They are simple, robust and particularly appropriate for the cell center finite volume schemes. In 
the following, the continuity criterion is first introduced in an heuristic way and then proved via the first principle of 
non-reflecting and plane wave propagation [1], Their limitations are also briefly discussed. 

For simplicity, the flow is assumed to be continuous near the boundary, i.e., with no shock or contact discontinuity. 
Later, it will be shown, as a discontinuous wave (e.g. a shock) may be decomposed by Fourier integral or Fourier 
series into a series of plane waves of different wave numbers, the NRBC may still be valid for it. 
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Figure 1: The NRBC criterion in E 3 . 

In the following context, U = (111,112, 113,114, uf) T = ( p , pu, pv, pw, pe) T , V = {vi,V2,v 3 , V 4 ,v$) T = ( p , u, v,w,p) 

and W = (w 1, 'W 2 • W 3 , W 4 , w 3 ) 7 are employed to represent respectively the conservative flow variables, the primitive 
flow variables and the characteristic variables. 

2.1 A heuristic approach 

We begin with the behavior of the characteristic variables W = [w\,W2, w 3 , w 4, w 3 ) T across an artificial boundary 
surface element As as time elapses. As also represents the interface between a boundary cell and its corresponding 
ghost cell. Note that any spatial boundary is a cylindrical hyper-surface in the space-time E 4 space (Fig. 1). 

There are various techniques to treat the NRBC based on characteristics theory. For instance, a well-known 1- 
D flow NRBC treatment by Enquist and Majda [3], and Hedstrom [7] is the requirement that the local perturbation 
(disturbance) along incoming characteristics be made vanish at the boundary (see [8], p.370). Let W = (w\, w 3 , w 3 ) T 
be the 1-D characteristic flow variables, the above requirement states that for those k such that the corresponding 
characteristic enters the computational domain through the artificial boundary: 

A w k = 0 (3) 


Let Pi and P e be respectively an interior point and an exterior point of the computational domain in the E,\ space. 
Both Pi and P e lie in the neighborhood of O, a point on As (Fig. 1). Then, 

Awk = Wk(Pi) — w k (Pe) = 0, for selected k. 

When Pi tends to O from the interior of the domain and P e tends to O from the exterior of the domain, the NRBC Eq. 
(3) becomes: 

lim w k (Pi) = lim w k (P e ) = w k (0) (4) 

Pi^O Pe — y O 
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for those selected ks. Therefore, the usual NRBC treatment is formally interpreted as a continuity problem of w k 
across the boundary surface. 

Secondly, for the NRBC of multi-dimensional flow, we formally extend the continuity of w k across the boundary 
surface to all ks rather than just a selected number of ks: 

lim Wk(Pi) = lim Wk(P e ), for all k. (5) 

Pi^O Pe — >0 

From the local equivalence of the characteristic variables W and V or U (see e.g. Hirsch[8], p. 155-156), Eq. (5) is 
equivalent to the continuity of V or U across the artificial boundary surface: 


^lim^ Uk(Pi) = p li m o Uk(P e ), for all k (6) 

or 

lim v k (Pi) = lim v k (P e ), for all k (7) 

Pi^O Pe — >0 

It is advantageous to switch from the continuity relation (5) of W to that of V or U, (6) or (7), since the latter can be 
treated in an easy way. At this stage, the continuity criterion of NRBC is heuristically inferred. 

2.2 the first principle of non-reflecting 

The above heuristic approach shows intuitively how the continuity criterion can be inferred from the characteristics 
NRBC. In the following, the continuity criterion will be proved, but first of all, a definition of ‘non-reflecting’ is 
needed: 

Definition: An artificial boundary is said to be non-reflecting if the solutions of the hyperbolic p.d.e.s (partial 
differential equations) in the domain interior and domain exterior are mutually solution continuations of each other 
across the boundary surface. 

In other words, both these solutions are but part of the solution in the entire domain. Hence, any point at the 
artificial boundary may be considered as an interior point in the entire domain and there is no reflection. This is the 
first principle of non-reflecting adopted in the present paper. If, in the particular case that simple wave-like solutions 
exist: 

V = Ve w , 

where i = \/— 1, then both the ‘amplitude’ V and the phase 6 are required to be identical across the artificial boundary. 


2.3 Plane wave solutions of the Euler/N-S equations 

The plane wave solutions are based on the Cauchy’s method of Fourier Integral (see Courant and Hilbert[l], pp.210- 
21 1). Consider the Euler/N-S equations in non-conservation form: 


dv ~dv ~dv ~dv ~ 

-dt +A &f + B ^ + c ^~ Q ’ 


( 8 ) 
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where A, B, and C are the jacobian matrices and functions of V. and the viscous terms may be included in the source 
term Q. 

In the neighborhood of a point O (xo, to) at the artificial boundary, ( 8 ) can be locally linearized by setting the 
jacobians A, B , C to their values at (). Assume a plane wave with the form 

V = Ve i6 , (9) 


and substitute in ( 8 ). Here 6 = 0(x, t) = x • k — ojt is the phase of the plane wave, i = \/—l. 9 = const, stands for 
a characteristic surface or wave front [ 1 ], ( 8 ) then becomes: 


where / is the 5x5 identity matrix. For any given wave number vector k = ( k x , k y , k z ), from the hyperbolicity of (8), 
real eigenvalues ojs exist such that | K — oil \ - 0 (dispersion relation), and so also 5 linearly independent eigenvectors 
r m , m=l,2,...,5 such that ( K — (jj m I)r m = 0, where the matrix K = k x A + k y B + k z C. 

As the eigenvectors r m are linearly independent, they may be used as the new basis for the 5 dimensional space 
and so V is converted into V = M - 1 V, where M -1 is the jacobian matrix of the coordinate transformation. Let 
V = v m r m , m= 1,2, 3,4,5 in turn, or, say, V = V m = (0, v m , 0, 0, 0) T in turn, where v m is the scalar function to be 
solved for, also let Q = q m Xm = Qm in turn. Then for ui = uj rn . ( 8 ) becomes 


9V ,, i - ,0V 
_ + 


r)W 

= Q e -* (x -k -u m t)_ 

ay oz 


Note that for each V m or Q m , only one component is non-zero (v rn or q m ). and the above equation is equivalent to a 
single scalar equation for v rn . Totally, there are five such equations. Also notice that changes with rn. Finally, the 
‘amplitude’ V may be expressed as 

5 

m= 1 

The same result was also given in [ 8 ] (Vol.l, p.153). Thus, ( 8 ) supports plane wave solutions in the form of (9) Thus, 
( 8 ) supports plane wave solutions in the form of (9). 

For more general waves other than the simple plane waves, as ( 8 ) is locally linearized in the neighborhood of the 
boundary point 0(xo, hi), they may be decomposed by Fourier series or Fourier integral with respect to wave number 
k and replaced by the superposition of plane waves. 

2.4 Proof of the continuity criterion 

With the presence of the artificial boundary s (hyper-surface), time-space E 4 is bisected into two portions, domain 
interior Dt and domain exterior D e (Fig. 2). Within each portion, the flow is governed by the same Eq. (1). From the 
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Figure 2: Sketch of the continuity criterion in 1-D flow, only a real component of V is shown. 


first principle of non-reflecting and plane wave propagation, it can be shown that locally there is no reflection at the 
artificial boundary surface s if the continuity criterion Eq. (7) (or Eq. (6), or Eq. (5) ) is satisfied. 

Proof: Consider the non-conservation form (8): 


dV ~dV ~dV 

~dT + A ^ + B ^ 


~ dV 

Cjr- = Q 
oz 


As a result of the continuity condition (7), an admissible given set of V at the boundary s may be used as a common 
boundary condition to solve separately for V,; and V e in their corresponding subdomains 1), and D e , which are 
separated by the artificial boundary s. (Note that generally, the admissible V given at s should be identical to the 
solution of V over the entire domain, see Appendix ). Here, Vj and V e are respectively the solutions of (8) in D-, 
and D e . Let V be the solution of (8) over the entire domain. Due to the uniqueness of solution for well-posed initial- 
boundary value problems ( see Appendix ), V,; is identical to V in Di and V e is identical to V in D e . Therefore, in a 
neighborhood of s, Vj and V e are mutually a continuation of each other across the boundary s and hence there is no 
reflection according to the definition in §2.2. 

To be more specific in terms of plane wave propagation, from §2.3, (8) supports plane wave solutions. As a general 
wave solution can be locally written as a superposition of the plane wave solutions by Fourier integral, it suffices to 
consider only the behavior of a single plane wave solution in the form of (9) at the artificial boundary s. 

Let O (x 0 , t 0 ) be any point at the artificial boundary s, then (8) can be locally linearized in the neighborhood of O, 
i.e.,/1, B and C are frozen at V Q as in §2.3. For any given wave number vector k, from the continuity criterion (7), 
V Q = Vj = V e , all A, B and C remain the same across the boundary and so also the eigenvalues cu* = w e = ui. At 

0, the plane waves Vj and V e share the same k,x, oj and t, and hence the same phase 6. Again, from the continuity 
criterion (7), Vj = V e , or 

y ie i( k » x -^i) _ V e e !; ( k,x ~“ t ), 

1. e., Vj = V e . Thus, the plane waves Vj and V e share the same amplitudes too. Therefore, Vj is completely identical 
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boundary 



Figure 3: Numerical treatment of NRBC in £ 3 . 

to V e in a neighborhood of O at the boundary (in terms of phase and amplitude), and there is no reflection at O across 
the artificial boundary surface. 

The continuity of V, or U or W (Eq.(5)- Eq.(7)) across the boundary surface is thus the basic criterion of NRBC 
adopted in the present paper. In §3.2, the numerical NRBC (Type I) is constructed based on the continuity criterion. In 
§4, the relation between an NRBC and the flux balance across the boundary surface is established. Such relation leads 
to another absorbing NRBC (Type II). 

3 The numerical treatment of the continuity criterion 

Fig. 3 illustrates a 2-D (in £ 3 ) NRBC treatment. Let AABC be a boundary cell centered at P, with the side 
BC coincident with the artificial boundary surface. A BCD is the ghost cell centered at Q, sharing the boundary 
edge BC with AABC. Let O be the centroid of the boundary surface element BCC'B' . The limiting process 
of lim p j _ > o U(Pj) is equivalent to extrapolating U from the interior node P to O by Taylor expansion. Similarly, 
limp e _>o U(P e ) is equivalent to extrapolating U from the exterior ghost node Q to 0 by Taylor expansion. 

Although theoretically, (7) implies up to C°° continuity, in numerical approximation, only low order continuity such 
as C ° , C 1 or C 2 , etc. can be achieved. Since a plane wave solution (9) is based on parameters such as its amplitude 
and phase, the numerical approximation of V is required to be at least C 1 continuous at the artificial boundary in order 
to be consistent with the physical solution. Taking a 1-D version of (9) for example, the C 1 continuity requirement is 
explained as follows. 

It suffices to consider only a scalar component of V, say, the first one p. After discretization, the (artificial ) 
boundary surface element center O is used to represent the entire surface element As. Then, from the continuity 
criterion, approximately, it can be inferred that the ‘amplitude’ (p 0 )i = (p 0 )e = Po , where the subscripts o, i and e 
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stand respectively for the surface center O , domain interior and exterior. Approximately, at the boundary surface As, 


Pi = p 0 S k ' X °~ U ' t °') =Pe = Po 


0 i(k e X o —UJe'to) 


( 10 ) 


Note that numerically the C° continuity result ( 1 0) provides no information about the wave number k and the frequency 
u>. With the presence of phase error, numerical reflection may still occur. However, if the numerical continuity is 
enhanced from C° to C 1 , i.e. 


( Pi)x — t kjpj — (Pe)x — ikePei ( Pi)t — i^iPi — ( Pe)t. — tCU e /? e , 

then ki = k e and u>i = u e and there is no phase error. 

In constructing the NRBCs, although any one of U, V, and W can be used, U is selected in the present paper. 
Therefore, in addition to U p , the space and time gradients of U at P, namely, U x , U y , U z and Ut are also required. 
The resulting linear Taylor expansion (C 1 continuity) yields better accuracy and is consistent with the NRBC criterion. 
The NRBC at the ghost node Q now turns out to be a problem of how to define U and its gradients at Q so that the 
flow is C l continuous at the boundary surface (represented by ()). 

3.1 Examples of NRBC - the extrapolation technique 

For the Type I (outflow) NRBC, under a mirror image assumption explained later, it is found that a simple extrapolation 
technique works well. 

First, an example of NRBC in P 3 (2-D space) for triangular mesh is illustrated. As shown in Fig. 3, assume A ABC 
is a triangular boundary cell with the edge BC lying on the boundary and conveniently parallel to the (/-axis. Define a 
ghost node D as the mirror image of the triangle vertex A with respect to BC. Then A ABC and A DBC are mutually 
mirror images of each other (the mirror image assumption). At time step n, conservative variables U are given at the 
cell center P of A ABC. Then, the NRBC (labeled as Type I) at the geometrical center Q of the ghost cell may be 
defined as: 

(U) Q = (U ) P , (U X ) Q = (U x )p = 0, (U y ) Q = (U y )p. (11) 

Apply linear Taylor expansions to domain interior and exterior separately: 

(U o) interior = Up + (jJO 2/p)(Uy)p T" l/2At(Ut)p, (U o) exterior ~ Uq T ( I/O 3/Q)(Uy)Q H - l/2At(Ut)Q 

hence, (U G ) exterior — (Uo) interior • 

Here, the subscripts P, Q and O of x and y denote the corresponding coordinates of P, Q and O , and from (1), the 
time partial derivatives (Ut)g = (Ut)p can be directly obtained. Thus, U is C l continuous at O across the boundary 
surface element, ( 6 ) is satisfied and the boundary surface element is non-reflective. 
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Figure 4: Control volumes (CVs) in for compact updating. Left: boundary cell A ABC and the corresponding 
hexagon CV base ASBQCRA; right: quadrilateral boundary cell ABCD and its corresponding octagon C V base 
AT BQC RD S A; li. S. T are centers of neighboring cells and BC the boundary. In any case, quadrilateral PBQC is 
a portion of the CV base. 

In a consistent way, for 3-D flows, under the same mirror image assumption on the ghost cells, the following 
extrapolations are valid NRBCs with C 1 continuity: 

U Q = Up, (U x )q = (U x )p = 0, (Uy) Q = (Uy)p,(U Z ) Q = (U z )p, (12) 

or 

U Q = U p + Ax(U x )p, (U X ) Q = (U X ) P , (U y ) q = (U y )p, (U Z ) Q = (U Z ) P , (13) 

where Ax = xq — xp. 

As demonstrated in the examples in §6 and §7, this Type I NRBC works well for either supersonic or subsonic flows 
at the outflow boundary, but it should be noted that: 

1. (12) or (13) is but apossible selection under the mirror image assumptions, there are otherforms of NRBCs based 
on (5) - (7); 

2. The extrapolation technique utilizes the nearby U P data to approximate the admissible U data at the artificial 
boundary, which is not an unreasonable choice, but there is a danger that the solution could drift away from the 
true solution (see Appendix ) after long time steps marching. A remedy is to incorporate the Type I NRBC with 
other physical boundary conditions (e.g. back pressure, etc.) 

The extrapolation technique is not a new idea, there are many successful examples, such as [9] and [10]. In [10], it 
was shown in a different way, that the extrapolation NRBC also works well even for subsonic flows. 

3.2 Numerical procedures of NRBC 

The implementation of the NRBC is incorporated in the numerical procedure and may be summarized in the following 
steps: 
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control volume 



Figure 5: Fluxes balance on an internal face Si in E 3 . 

(i) Based on the flow data at the boundary cell center P , i.e. Up and its spatial gradients (slopes) U x , U y and U z , 
determine the flow data Uq as well as its gradients U x , U y and U z , i.e. the NRBC at Q, as described in §3.1. 

(ii) Update U at boundary cell center P to the new time level by the conservation laws (2). In order that U be C 1 
continuous across the artificial boundary, the updating procedure must be carefully designed to take account of the 
accuracy of surface flux calculation. Here a compact updating procedure described in §3.3 is recommended. 

(iii) After U at all the interior cell centers of the computational domain are updated, evaluate the new spatial 
gradients U x , U y , U z at the boundary cell center P by finite difference. Certain flux limiters might be needed in the 
process. For multi-dimensional flows, a linear equation system is required to solve for the gradients. 

(iv) Repeat steps (i) - (iii) and march in time. 

3.3 Compact updating - a viable time marching technique 

The updating procedure recommended here is identical to the one used in the recent CE/SE method (Chang et al 
[12,13]) and similar to the NT (Nessyahu-Tadmor ) scheme [14] in 1-D flow. 

The purpose for compact updating is to achieve high accuracy ( C 1 continuity) from a small cell stencil (e.g. the 
stencil formed by the immediate neighboring cells). Therefore, not only U but also its gradients U x , XJ y , U~ and U* 
are required. The compact updating is capable of maintaining C 1 continuity of U across the artificial boundary. 

Consider a triangular cell A ABC in E 3 (Fig.4, shaded area). Let P and Q, R, S be respectively the cell centers of 
A ABC and its neighbouring triangular cells. Most of the finite volume schemes use the space-time cylinder A ABC 
as the control volume (CV) for updating V or U at P. Surface flux along, say, AB is obtained by extrapolation from 
S across the cell to the surface center along AB. Another treatment is to include S in the CV and replace the flux 
along AB by the fluxes along A S and SB. Since S sits right on both surfaces along .4 S and SB, no extrapolation 
across the interior of the CV is involved. For the triangular cell A ABC in Fig. 4 the CV turns out to be a space-time 
hexagon cylinder in E 3 based on ASBQCRA. By applying the integral conservation laws (2) to the CV, updating 
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flow data at P based on the flow data at a compact node stencil of Q , R, S is now completed. Fig. 4 also demonstrates 
that for quadrilateral mesh cells the CV is an octagon cylinder in E$ (2-D space). 

The surface fluxes for the CV surfaces passing through a vertex, say, cell centers Q, R, S, can be evaluated by first 
extrapolating U along the surface to their corresponding surface centers by linear Taylor expansion, calculating flux 
functions F, G and H, and then incorporating the surface unit normal vector and computing the fluxes. For high space 
dimensions, the CV s are geometrically more complicated. More details including the updating of U x , U :7 , and U z 
can be found in [12,13], 

It should be noted that the compact updating is suggested for boundary cells only. For interior cells, one may resume 
to use other finite volume schemes. 


4 The absorbing NRBC and the matched layer 

Based on the same continuity criterion, another NRBC (Type II) is introduced in this Section. With this NRBC, an 
interesting extra advantage is, a thin matching layer between the interior flow and the prescribed boundary condition 
will be formed automatically should any discrepancy is developed between them. 


4.1 Relation between an NRBC and the flux balance across the boundary surface 

In this subsection, the relation between two statements is established. The first one states that the incoming fluxes at 
the artificial boundary surface are equal to the outgoing fluxes, or fluxes are balanced across the boundary surface. 
Here, the outgoing flux is defined as a portion of the left hand side of (2): 


ds. 


to = 1 , 2 , 3 , 4 , 5 , 


Jas 

where AS* denotes the artificial boundary surface (element). The second statement is that the boundary is non- 
reflective. Such a relation leads to the Type II NRBC. Numerical implementation of the Type II NRBC is straightfor- 
ward and will be illustrated in §4.2. 

Let V be any control volume in the P,\ space intersected and divided into two portions V\ and Vi by an internal 
surface Si. Let E-| , S 2 and o\, <ji be the fluxes around the surfaces of Vi and Vi respectively. Here o\ and cr 2 are the 
outgoing fluxes at the interface S, for V\ and Vi respectively (Fig. 5 ). Then, the following lemma holds: 

Lemma 1: For a control volume V in the E\ space, fluxes passing through any of its internal surface Si are 


balanced, i.e. cri + <r 2 = 0. 


Proof: Apply (2) to V, V\ and Vi separately. We have 


Si + S2 — / QmdV, Si T <j\ / Q m dV§ S2 + (t 2 — / Q m dl'. 

Jv J \ 1 Jv 2 

Add the second and the third equations and then subtract the first one, we have a\ + 01 = 0. 

In other words, across any interior surfaces, no extra ‘source’ should be generated. 
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artificial boundary surface 
element ds 



Figure 6 : Fluxes balance across a boundary surface element in a control volume (the above figure) in E 3 . 


Next, we show that locally the continuity of flow variables (Eq. (5), ( 6 ) or (7)) can be inferred from flux balance. 
Hence, flux balance suffices for an NRBC. 

Consider an element ds of the cylindrical spatial boundary surface in E 4 . As shown in the control volume in Fig. 6 
(in E 3 ), ds is centered at O. Assume the outgoing unit normal vector at O is ni = (n x ,n y , n z , 0) T , then, the incoming 
unit normal is ri 2 = {—n x , —n y , —n z , 0) T . The outgoing flux <j\ for the inner boundary surface: 
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Let the outgoing boundary surface flux vector be: 


ds[n x F + n y G + n z Il\. 


L — [rt^F T Tiy G T rirH], 


then 


L = ai/ds. 


Let V = (p. u, v, w, p) r be the primitive flow variable vector. L may be considered as a non-linear vector function of 
V. The jacobian matrix 4^ has the eigen values (e.g. see Hirsch [8], p.177): 


Ai = A 2 = A 3 = u ■ n x + v ■ n y + w ■ n z , A 4 = Ai — c, A 5 = Ai + c. 

where c is the speed of sound. If none of the eigenvalues vanishes at O, the jacobian ^ is non-singular, and an inverse 
vector function, the primitive flow variables V as a vector function of the surface flux L exists. Thus, 


NASA/CR— 2003-212387 


13 



boundary surface 



Figure 7: NRBC at inflow in E 3 . Inflow data are given at the ghost cell center. 

Lemma 2: If the jacobian is non-singular, then, locally, the primitive variables V are uniquely defined by the 
flux vector L at the center O of the boundary surface element ds. 

Proof: Assume L(Vi) = Li and there is another V 2 in the neighborhood of Vi, such that L(V 2 ) = Li. By a 
linear Taylor expansion, 

AT 

L(V 2 ) - L(Vi) = (— )v= Vl (V 2 - Vi) + 0(1 V 2 - Vij 2 ) = 0 
aV 

Hence V 2 — Vi = 0orV2=Vi since the first order term cannot cancel with the second order term. 

From Lemma 2 and the NRBC continuity criterion (7), it is inferred that locally, under the condition that the jacobian 
is non-singular, the following lemma holds: 

Lemma 3: For hyperbolic conservation laws of gas dynamics, an element of the artificial boundary surface is non- 
reflective if its outgoing fluxes and incoming fluxes are equal (balanced). Lemma 3 states that conditionally the flux 
balance across a boundary surface element is a sufficient condition for NRBC. 

With the concept of flux balance, there is also an intuitive interpretation of the commonly used terms for NRBC 
such as ‘transparent’ and ‘absorbing’. If the flux 100% passes through the boundary surface element (or the flux is 
balanced), the boundary surface is said to be ‘transparent’ or ‘absorbing’ to the fluxes. Chang et al [11] were the first 
attempting to explain the 1-D NRBC using flux concepts. 

4.2 Implementation of absorbing boundary condition 

The lemmas in §4.1 can be easily applied to construct the Type II NRBC. In case that at the ghost cell center nodes 
flow variables must be specified as the given values (e.g. at the inflow boundary, Fig. 7), another type (Type II) of 
NRBC - absorbing NRBC arises. A control volume V across the boundary surface is needed to apply the divergence 
theorem (2). As shown in Fig. 7, the ghost cell center lies outside of the domain, and the boundary surface element As 
is an internal surface of the control volume V, then, from Lemma 1 and Lemma 3, As is automatically a non-reflective 
boundary surface element. 
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The Type II NRBC states that no extra condition of NRBC is needed with the prescribed flow boundary conditions at 
the ghost cell center. This Type II NRBC is flexible and valid with practically any cell shapes or configurations. Unlike 
the Type I extrapolation NRBC, ghost cells are not required to be mirror images of the corresponding boundary cells. 
However, it is noted that U is still required to be C 1 continuous across the artificial boundary surface. Consequently, 
the compact updating is still recommended for the Type II NRBC. 

In many cases, over a long time period, the (subsonic) flow may develop a flow near the inflow boundary that 
is different from the imposed boundary conditions, or the boundary condition becomes ‘overdetermined’ . Then, a 
‘matched layer’ may be formed between the computed interior flow and the imposed flow boundary condition (see 
Fig. 18). The situation is somewhat similar to that in the PML (perfectly matched layer) method [5,6] except that the 
same set of governing equations are still valid in the layer. The shock-capturing numerical scheme should be able to 
quickly resolve this non-physical discontinuity in a few cells. 

5 Discussions on the NRBCs 

In practical applications, due to discretization, approximation and lack of information in the domain exterior, there 
are limitations for both Type I and Type II NRBCs. 

For Type I NRBC, as mentioned in §3.1, extrapolation technique might lead to drifting or deviation from the true 
solution because no flow data outside the outflow boundary is available. In addition, under the mirror image assumption 
on the ghost cells in §3.1 and the assumption that the boundary surface As is normal to the x axis, the NRBC (12) 
implies that the NRBC continuity criterion (5 - 7) are satisfied at any point on As. However, as explained in the 
following, due to discretization and the possible consequent phase error, the accuracy of the NRBC could be degraded. 

Consider the Fourier mode in a plane wave solution (9), i.e., e *( k,x - wt ). Here, 9(x, t) = k • x — cut is the phase 
of the wave mode, with k being the wave number vector in the propagation direction. Generally, the direction of k 
may or may not be the same as the flow direction. 0(x, t) = const, stands for a wavefront ( or a characteristic surface, 
see e.g. Hirsch [8], p.150, Courant and Hilbert [1] ). After discretization, the center O of As (Fig. 3) is employed 
to represent the entire As. Then how much phase error is introduced to the Fourier mode by the discretization? Let 
x = (x, y, z) be the position vector of any point on As and xo the position vector of the center (). For clarity, assume 
time t is held unchanged. Then, after discretization, the phase error A 0 due to replacing x by xo is: 

A 9 = k • x — k • xo = k • (x — xo) (14) 

note that Ax = x — xo lies on As, A 9 = 0 when k is normal to As. Therefore, for Type I NRBC, the best result is 
obtained when the wave propagation direction is normal or only slightly oblique to the boundary surface. Otherwise, a 
phase error of order O(Ax) may be introduced. It deteriorates the accuracy of NRBC and causes numerical reflection. 
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Figure 8: Comparison of exact and numerical results at the domain boundary for the 1-D Gaussian pulse problem. 



Figure 9: Propagation of 1-D Gaussian pulse at the domain boundary. 

For Type II NRBC, in addition to the similar phase error of Type I NRBC, there are other restrictions too. In §4.1, 
Lemma 3 is conditionally valid because it is based on the one to one correspondence between the boundary surface 
flux vector L and the primitive flow variable vector V. The latter relies on the non-singularity of the jacobian In 
case that ^ is singular, Lemma 3 may fail. In addition, Lemma 3 is valid only locally. Globally, the relation between 
L and V involves a quadratic equation, the vector function V (L) could be multi-valued. This could break the one to 
one correspondence between L and V globally and lead to the failure of Lemma 3 in the global sense. 

In a nutshell, in reality, there are various situations that the Type I or II NRBCs can only be partially or approximately 
implemented, causing spurious reflections at the (artificial) boundary or deviation from the true solution. In practice, 
an effective remedy is to impose a buffer/sponge zone between the boundary and the interior domain. Although the 
same governing equations (1) or (2) are employed in the sponge zone, numerical damping is highly increased to 
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diminish the wave/disturbance amplitude before it reaches the boundary and to minimize the spurious reflection. An 
example is depicted in §8. 


6 Numerical examples for outflow NRBC 

In this Section and §7, the effectiveness of these NRBCs is demonstrated in numerical examples in one, two and 
three dimensional spaces. 

In principle, any finite volume scheme can be used with the above NRBC if it can be manipulated at certain high 
accuracy. Here the recently developed space-time conservation element and solution element (CE/SE) method [12,13] 
is chosen for computing the examples since the compact updating is a standard procedure in the scheme, making 
application of the NRBCs straightforward and effective. Full details of the method are described in [12,13], The Type 
I NRBCs used with the CE/SE method are identical to (12), with possible minor modification according to the grid 
layout. 


6.1 Propagation of a 1-D Gaussian pulse 

Consider a scalar initial value problem: 


du du 
~dt + ^c = 


over the range —20 ^ x ^ 450, with a Gaussian pulse u = 0.5 exp [— (7n2)(|) 2 ] at t = 0. This is one of the 
benchmark problems of the 1st CAA Workshop [15]( Category 1, Problem 1). The exact solution given there is: 

,x — t . 


u = 0.5 exp 


-{In2){—Y 


In this example. Ax = 1 is chosen and At = 1 is based on CFL number = 1. With CFL number = 1, and other 
parameters e = 0 and a = 0, the CE/SE scheme yields a numerical result which is identical to the exact solution in 
the interior of the domain —20 ^ x ^ 451. Thus, performance of the outflow Type I NRBC at x = 451 can be easily 
validated. At t = 450, the Gaussian pulse is passing through the outflow boundary x = 451, where the Type I NRBC 
(12) with appropriate modification to 1-D flow is imposed. The table in Fig. 8 lists the exact solution and the CE/SE 
result from x = 430 to x = 451. The result is also plotted in Fig. 9. It is seen that they are completely identical (up to 
10 decimal places) and there is absolutely no reflection, although the grid is rather crude with Ax = 1. 


6.2 2-D free shear layer instability and vortex roll-up 

The problem considered here is identical to the inviscid free shear layer instability problem considered in [16] (Fig. 10). 
The background mean flow consists of a fast stream (supersonic) in the upper half domain and a slow stream (subsonic) 
in the lower half domain. The two parallel streams are connected by a continuously changing shear layer of the 
hyperbolic tangent profile. 
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Figure 10: Free shear layer instability problem, u\ = U\ = 1, v\ = 0, pi = 1/3.15, p\ = 1, M\ = 1.5, U 2 = 
U 2 = .7391304, V 2 = 0, P 2 = 1/3.15, P 2 = 0.5405405, with subscripts 1,2 denoting the fast and slow streams 
respectively. 



7 - V 






200x20 domain, 200x100 grid 
(same contour levels) 

— 100x20 domain, 100x100 grid 


comparison of numerical results with different computational domains, 
showing effectiveness of outflow NRBC 



Figure 11: Contours for long and short domains, showing effectiveness of the outflow NRBC. 
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p along x=y,ep=.1 ,A1=.6,at t=100(Cat.4.1) 



Figure 12: An acoustic pulse above a solid surface passing through the outflow boundary. 

In the test, two computational domains are chosen. The first one is 0 ^ x ^ 200 and — 1 0 ^ y 10, with a grid of 
200 x 100 uniform cells. The second one is 0 ^ x ^ 100 and — 10 ^ y ^ 10, with a grid of 100 x 100 uniform cells. 
Both cases have exactly the same grid cell sizes, time step size At = 0.15, and parameters e = 0.2, a = 0. They are 
both run for 4000 time steps when the spatial instability is fully developed. To ensure that the instability waves and 
vortex roll-up develop quickly, a large perturbation amplitude of 0.02 at the most unstable frequency is chosen for the 
eigenfunctions. At the outflow boundaries, the Type I outflow NRBC is used. Figure 1 1 demonstrates snapshots of 
the isobars and isopycnics in the two cases. These contours are observed to be almost identical to each other in their 
common domain portion. The contours in the short domain seem as if they were a piece chopped off from the longer 
one. This shows that the outflow NRBC in this case is nearly perfect. 

6.3 Acoustic Pulse Propagation 

This problem is a typical subsonic wave propagation problem [15]. The computational domain in the x-y plane is a 
square with —100 ^ x ^ 100, and 0 ^ y ^ 200. A uniform 201 x 201 (triangulated) grid is used with Ax = Ay = 1 
. Initially, a Gaussian acoustic pulse is located at the lower portion of the domain (x = 0, y = 25), with a mean flow 
of Mach number M = 0.5 in x direction and a solid wall at the bottom. At the other three boundaries. Type II NRBC 
is used. By choosing a small amplitude factor 6 = 0.001, the Euler equations are practically linearized. Fig. 12 shows 
the isobars at t = 100 and the comparison between numerical and analytical results for density along the line x = y 
for 0 ^ x ^ 100. Although the wave propagation direction is oblique to the outflow boundary, only tiny reflection is 
observed from Fig. 12. Another 2-D example with completely subsonic flow will be illustrated in §7. 
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Figure 13: Sketch of the rectangular jet, aspect ratio 5, jet Mach number Mj = 1.6, L = 1 6 1) . IF = 1613 and 
H = 5.6 D. 

6.4 3-D rectangular jet flow 

Fig. 13 is the sketch of an underexpanded rectangular jet in 3-D space. The rectangular nozzle protrudes into the 
computational domain by l = 2D, D being the width of the jet. The unstructured mesh consists of about 1.7 million 
tetrahedral cells. At the inlet plane, ambient (stationary) condition is specified. Jet flow at higher pressure is specified 
at the nozzle exit. All the rest boundaries are either Type I or Type II NRBC. Fig.s 14 shows snapshots of the isobars 
and v velocity contours on the cross sectional mid-planes after running 60,000 time steps. Fig. 15 demonstrates the 
3-D pressure iso-surfaces. No visible reflection is observed. 

6.5 Influence of the NRBC to the numerical accuracy 

In §3.1 and §5, it is commented that since the admissible V at the artificial boundary is not available and has to be 
approximated by extrapolation, this will eventually lead to deviation or drifting of V from the true solution. But in 
realty, in many cases, the influence from the boundary to the interior flow is small. An example is shown here. 

Fig. 16 demonstrates the sound intensity level computed based on two domains with different lengths but same 
width for the 2-D Mach radiation problem in the 3rd CAA Workshop (Category 5) [18]. The short domain has a grid 
of 289 x 144 nodes. The only difference is that the longer domain has 30 more uniformly distributed nodes added 
in the x direction downstream. Fig. 16 shows the sound intensity levels ( square of r.m.s. p' - pressure fluctuation ) 
along the line y = 10 at t = 400 (40000 time steps or 28 periods ) for both domains. It is observed that the maximum 
difference is about 2 x 10~ 10 , far below the discretization error, thus is negligible. This case also demonstrates the 
relative drifting side effect of the extrapolation Type I NRBC, but the error is acceptable. 


NASA/CR— 2003-212387 


20 




v- contours on the mid plane (narrow side) 



v-contours on the mid-plane (wide side) 


Figure 14: v velocity contours on the mid-planes with mesh background. No visible reflection is observed. At the 
outflow boundary, flow is supersonic in the jet core and then becomes subsonic across the thick shear layer. 


NRBC 



Figure 15: Pressure iso-surfaces, no reflection is observed. 
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Figure 16: Comparison of sound intensity levels from computational domains of different lengths, showing that the 
outflow NRBC has negligible influence to the interior. 

7 Numerical example for the absorbing boundary condition 

As stated in §4.2, for an absorbing boundary condition such as the inlet, the prescribed inlet conditions are already 
good enough for an NRBC. Fig. 17 illustrates the instantaneous isobars for a cavity noise problem. 42, 000 triangulated 
structured cells are used. The problem is a M=0.8 flow past a rectangular cavity of aspect ratio of 6.5. At the cavity 
walls, no slip boundary condition is imposed. Due to vortex shedding and acoustic feedback at the cavity edges, strong 
nonlinear acoustic waves are generated and propagate in both upstream and downstream directions [17]. Fig. 18 is an 
enlargement of Fig. 17 around the inlet area. The details of the contours at the matched layer is revealed. It is observed 
that there is no spurious reflection and that the matched layer is about 4-5 cells’ thick. The matched layer in the Type II 
NRBC is somewhat similar to that of the PML (perfectly matched layer) method [5,6] in that the difference diminishes 
quickly within this layer. But the layer arises automatically and there is no need to solve a new set of equations in the 
layer or to impose any conditions other than the prescribed inflow physical conditions. 

For NRBC at the top of the computational domain, a Type I NRBC may be modified from (12) by simply exchanging 
the axes x and y: Uq = U p , (U y )Q = (U x )q = 0, where as before, P and Q are respectively the boundary cell 
center and the ghost cell center at the top boundary. From Fig.s 17-18, even when the acoustic wave is oblique to 
the boundary, there is still no visible reflection. A Type I NRBC is applied to the outflow boundary, still, no visible 
reflection is found for the subsonic flow. 
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Figure 17: Isobars snapshot for a cavity noise problem (Mach number M = 0.8), showing inflow NRBC and its 
absorbing property. No visible reflection is found at the top and the outflow boundary 



Figure 18: Details of the contours at the inflow boundary, showing the matched layer at the inlet and its spreading over 
the grid. 
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Figure 19: Mach radiation from a M = 2 axisymmetric jet (without buffer zone), showing severe spurious numerical 
reflection. 



Figure 20: Mach radiation from a M = 2 axisymmetric jet (with buffer zone but not shown), showing a clean acoustic 
field. 


8 Application of the buffer/sponge zone 

Generally, the NRBCs amount to little reflection. However, there are situations that they may fail and discernible 
reflections occur, such as the cases discussed in §5. A simple but effective remedy is to add a buffer/sponge zone 
between the core domain and the boundary. In the buffer zone, the same governing equations are still used, except that 
numerical damping is highly increased. Here, the cell size in the buffer zone may grow rapidly (e.g. exponentially ) 
and create larger numerical damping. Typically, the number of cells in a buffer zone may vary from a few to 20. An 
example of a 2-D axisymmetric jet with Mach radiation from an externally stimulated shear layer is demonstrated. It 
is similar to the one described in the 3rd CAA workshop benchmark problems (Category 5) [18]. The domain size is 
3319 x 19 D with D being the jet nozzle diameter. 300 x 280 non-uniform rectangular cells are used before they are 
further triangulated. Initially, a Mach number M = 2 jet exists. At the center of the nozzle exit plane, a source is 
imposed and perturbs the jet flow with a small amplitude A = 0.001 at a Strouhal number St = 0.2. Mach radiation 
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is then triggered and gains its strength along the stream. At the outflow boundary, particularly at the shear layer, 
due to the staggered type mesh and the strong Mach waves with oblique angles, the Type I NRBC fails and spurious 
reflection is generated and propagates upstream (Fig. 19). After a buffer zone of 10 cells is added at the domain 
outflow boundary, the spurious reflection disappears and a clean mach wave radiation is shown (Fig. 20). The size of 
the buffer zone cells grows exponentially at a rate of 20% along the x axis. 

9 Concluding Remarks 

In the present paper, based on the first principle of non-reflecting, the propagation of plane waves and the hyper- 
bolicity of the Euler equation system are revisited, and then combined to derive the continuity criterion of NRBC. 
This continuity criterion reveals the nature of NRBC and it is possible to apply it to hyperbolic systems other than the 
Euler equations. Numerically, C 1 continuity is required at the artificial boundary to be consistent with the criterion 
and achieve non-reflection. 

Two basic types of NRBCs are described: the extrapolation (Type I) and the absorbing (Type II) NRBCs. Simple 
but effective C 1 continuity NRBCs are consequently developed. Generally, their performances are similar to those of 
the characteristics-based NRBCs. For shock-capturing schemes. Type II NRBC often forms a matched layer within 
a stencil of a few cells, should discrepancy develops between the interior flow and the boundary conditions. This 
matched layer is similar to the well-known PML (perfectly matched layer) in [5,6], Limitations of the NRBCs are 
also discussed. In particular, the Type I extrapolation NRBCs may cause solution drifting due to lack of information 
beyond the outflow boundary (see also [10]). A remedy is to incorporate the physical boundary conditions or to use a 
sponge zone (cf. Appendix). 

The diversity of various NRBCs in flow computations can never be overestimated. The purpose of the present paper 
is to show some guidelines in developing NRBCs that are simple and robust for practical computations. The compact 
updating procedure proves to work well with the NRBCs and provides C 1 accuracy in surface flux evaluation. But it 
is definitely not the only way to achieve non-reflecting effect. Different schemes may have different treatments, see for 
example, [9,10], Sometimes, a combination of the NRBC treatments may provide much improved results, such as the 
incorporation of the buffer zone. As a byproduct, we are now also able to explain why the NRBCs with the recently 
developed CE/SE scheme [12,13,16,17] are robust. 

The restriction set in §2 that the flow is continuous may be lifted, since on the one hand, for shock-capturing 
schemes, a discontinuity may be considered as a continuous wave with steeper gradient, and on the other hand, a sharp 
gradient may be decomposed by Fourier series or Fourier integral. But if the wave propagation direction is oblique to 
the (artificial) boundary, the NRBC may perform poorly due to reasons explained in §5. Chang et al [11] presented 
a 1-D example ( i.e. wave propagation direction normal to the boundary) showing how a shock passes through an 
artificial boundary without causing visible reflection. 
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Appendix: Discussions on the continuity criterion and the extrapolation technique 

Without loss of generality, consider a 1-D flow shown in Fig. 21. As a pure initial value Cauchy problem, V = 
V(x, 0) at t = 0 is given. This problem is well-posed and there exists a unique solution V = V 0 (x,t) over the 
entire domain D: -oo < x < oo, t ^ 0. If a boundary exists on the left inlet side, the problem becomes an 

initial-boundary-value problem, a boundary condition is required at the inlet boundary. 

Assume the artificial boundary locates at x = 0 and bisects the entire domain D into domain interior I), and domain 
exterior D e : D = Di + D e (see §2.4). In order that V,; and V e in §2.4 are respectively the unique well-posed solutions 
in subdomains Di and D e , an admissible common boundary condition they share at x 0 is: 

Vj(0,t) = V e (0,i) = V o (0,t). 

Here, V o (0, t) is the entire domain solution along the line x = 0; or equivalently, V o (0, t) may be obtained by the 
method of characteristics as sketched in Fig. 21 at point O at the boundary. 

However, for NRBC problems in reality, V 0 (,x, t). x Js 0; or V e (x. t) is totally missing (otherwise there is no 
need to investigate the NRBCs). In this situation, the extrapolation technique e.g. (12) is not an unreasonable choice 


t 

i 



Figure 2 1 : Continuity criterion in 1 -D flow. 
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for the NRBC continuity criterion. The initial-boundary value problem for V , is well-posed only when the flow is 
supersonic or is known to remain unchanged across the artificial boundary at x = 0 by a priori information. Otherwise, 
even though locally the extrapolation f 12) leads to non-reflecting effect, globally the ‘solution’ will keep drifting away 
and deviate from the true solution. The following are some remedies for practical numerical flow computations: 

[1] use a sponge (buffer) zone with highly increased numerical damping to filter away the wave ingredients in the 
flow; when the Type I extrapolation NRBC is applied to the outer outflow boundary, the flow is already uniform 
at the boundary. 

[2] incorporate the extrapolation with other physical boundary conditions (e.g. back pressure etc.). 

Notice that the present discussions do not apply to the Type II NRBCs. 
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